//version 0.4 (Nov 22, 2013).
// tempfile fixed from 0.1
// markout/touse fixed from 0.2
// precision (gen/egen double) fixed from 0.3
//TO DO: fix noadj option (Jul 2013).
cap program drop arsens

program define arsens,rclass
version 9


syntax varname [if] [in], Treat(varname numeric) Str(varname numeric) [Gam(real 1.0) Del(real 1.0) Verbose NOAdj Kverbose] 

if `c(changed)' == 1{
di "{error:Save data in memory before proceeding (note: clear option invalid for arsens)}"
exit 4
}

if `gam'<1{
di "{error:Enter a value of Gamma that is at least 1}"
exit 198
}

if `del'<1{
di "{error:Enter a value of Delta that is at least 1}"
exit 198
}

local mainfile= "`c(filename)'"

//begin quietly
quietly{

local rusemui=0
local rusesigmi=0

marksample touse
markout `touse' `treat' `str'

//for clarity in code:
local resp="`varlist'"

//b/c input is GAMMA, we need to covert to gamma;
//exp(gamma)=GAMMA--> ln(GAMMA)=gamma. 
//mutatis mutandis, delta


local udel=ln(`del')
local ugam=ln(`gam')





//stratum avg resp:
egen double stavresp=mean(`resp') if `touse'==1, by(`str')
//"aligned response":
gen double alresp=`resp'- stavresp if `touse'==1
//aligned rank:
egen double alrank=rank(alresp) if `touse'==1
//list stavresp alresp alrank treat resp strat



//adjusted  alrank

if "`noadj'"=="" {
quietly summ alrank if `touse'==1
replace alrank=alrank/(.5*`r(max)') if `touse'==1
}



//observed test stat (?)
egen double tv=total(alrank) if `treat'==1 & `touse'==1
egen double tv2=mode(tv)
local tsobs=tv2


drop tv tv2


//replaced-->gen __worb=`resp'

gen double __worb=alrank

//replaced-->gen double __rorb=`resp'

gen double __rorb=alrank


gen double __zorb=`treat' if `touse'==1

gen double __borb =`treat' if `touse'==1


gen double _gunobs=`ugam' if `touse'==1
gen double _dunobs=`udel' if `touse'==1


sort __worb
//note equivalent to "sort __rorb", and 
//probably irrelevant in any case.

gen double amuik=.

gen double muik=.
gen double sigmik=.

gen double  usemui=.
gen double  usesigmi=.
//ensure the above 4 are re-set as/if needed.
//(which they are b/c new file opened for every i-loop.)


tempfile tfile
save "`tfile'",replace

quietly summ `str' if `touse'==1
local strmax=r(max)

//START LOOP 1
forvalues i=1/`strmax'{ 
use "`tfile'",clear
marksample touse
keep  if `str'==`i' & `touse'==1

if `c(N)'==0{
continue
}

if "`verbose'"=="verbose" | "`kverbose'"=="kverbose" { 
noisily di "STRATUM: `i'" 
}



tab _gunobs 
local K_i=`r(N)'




//START LOOP 2
forvalues k=0/`K_i'{

if "`kverbose'"=="kverbose"{
noisily di "K: `k'"
}

capture drop arunmuik runmuik runsigmik

gen double arunmuik=0 in 1
gen double runmuik=0 in 1
gen double runsigmik=0 in 1


capture replace _gunobs=0 in `k'
capture replace _dunobs=0 in `k'
//replaced-->mata: arssub1 ("`resp'")
mata: arssub0 ("alrank")
mata: arssub0V ("alrank")
//return list
//tab _gunobs
} //END LOOP 2


local usemui=usemui
local rusemui=`usemui'+`rusemui'

local usesigmi=usesigmi
local rusesigmi=`usesigmi'+`rusesigmi'





} //END LOOP 1 


local dev=(`tsobs'-`rusemui')/(sqrt(`rusesigmi'))


local pval=1-normal(`dev')

use "`mainfile'",clear

ret scalar pval=`pval'
ret scalar dev=`dev'
ret scalar tsobs=`tsobs'
ret scalar gam=`gam'
ret scalar del=`del'
ret scalar summu_imax=`rusemui'
ret scalar sumsd_imax=`rusesigmi'


//end quietly:
}

di  _newline(1) "{result:{ul on}SIMULTANEOUS SENSITIVITY ANALYSIS FOR ALIGNED RANK TEST{ul off}}" _newline(2) ///
"{res:Gamma: `gam'}"  _newline(1) ///
"{res:Delta: `del'}"  _newline(1) ///
"{res:Expectation: `rusemui'}" _newline(1) ///
"{res:Variance: `rusesigmi'}" _newline(1) ///
"{res:Deviate: `dev'}" _newline(1) ///
"{res:max[p(t>=`tsobs'|u)]: `pval'}"




end

version 9
mata:
void arssub0 (transmorphic matrix input0a){
//input0a not used.


GUVEC=st_data(.,"_gunobs")
DUVEC=st_data(.,"_dunobs")
BVEC=st_data(.,"__borb")
WVEC=st_data(.,"__worb")
DEN1=0
//DEN1

infoB=cvpermutesetup(BVEC) 
while ((BVEC=cvpermute(infoB))!=J(0,1,.)){
//GUVEC
//BVEC
preDEN1=exp(quadcross(GUVEC,BVEC))
//preDEN1
DEN1=DEN1+preDEN1
//DEN1
}


DEN2=0
DEN2
infoW=cvpermutesetup(WVEC) 
while ((WVEC=cvpermute(infoW))!=J(0,1,.)){
//WVEC
preDEN2=exp(quadcross(DUVEC,WVEC))
DEN2=DEN2+preDEN2
//DEN2
}




arssub1 ("alrank",DEN1,DEN2)

}


end



version 9
mata:
void arssub1 (transmorphic matrix input1a, real vector DEN1, real vector DEN2){ 

RVEC=st_data(.,input1a)



infoR=cvpermutesetup(RVEC) 

 
while ((RVEC=cvpermute(infoR))!=J(0,1,.)){

st_store(.,"__rorb",RVEC)
arssub2 ("__zorb",DEN1,DEN2)

}


}


end




version 9
mata:
void arssub2 (string scalar input2a, real vector DEN1, real vector DEN2){ 


ZVEC=st_data(.,input2a)
RVEC=st_data(.,"__rorb")


//RVEC
//ZVEC

//RSUM=0

infoZ=cvpermutesetup(ZVEC)  
while ((ZVEC=cvpermute(infoZ))!=J(0,1,.)){
arssub3(ZVEC,RVEC,DEN1,DEN2)
//ZVEC
//RVEC
}
}


end



version 9
mata:
void arssub3 (real vector zV, real vector rV, real vector DEN1, real vector DEN2){
GUVEC=st_data(.,"_gunobs")
DUVEC=st_data(.,"_dunobs")
BVEC=st_data(.,"__borb")
WVEC=st_data(.,"__worb")

NUM1=exp(quadcross(GUVEC,zV))
NUM2=exp(quadcross(DUVEC,rV))

TESTST=quadcross(zV,rV)


AMUIK=(TESTST)*((NUM1)/(DEN1))*((NUM2)/(DEN2))
//AMUIK
//TESTST 
//NUM1
//DEN1
//NUM2
//DEN2

//note that this is the summand of MU_ik; not MU_ik itself.
//mu_ik is only calculated in arssub4.


//(note case sensitivity...)
st_store(1,"amuik",AMUIK)








stata ("arssub4 amuik")
}


end 




cap program drop arssub4

program define arssub4,rclass
version 9

syntax anything(name=numl)

gettoken one: numl




qui replace arunmuik=arunmuik + `one'
//note that this gives only one observation, which is what we want.
//(also one of runmuik, arunmuik is redundant, really).




end


version 9
mata:

void arssub0V (transmorphic matrix input0a){
//input0a not used.



//ensure __borb and __worb is still OK.
GUVEC=st_data(.,"_gunobs")
DUVEC=st_data(.,"_dunobs")
BVEC=st_data(.,"__borb")
WVEC=st_data(.,"__worb")
DEN1=0
infoB=cvpermutesetup(BVEC) 
while ((BVEC=cvpermute(infoB))!=J(0,1,.)){
//GUVEC
//BVEC
preDEN1=exp(quadcross(GUVEC,BVEC))
//preDEN1
DEN1=DEN1+preDEN1
}


DEN2=0
infoW=cvpermutesetup(WVEC) 
while ((WVEC=cvpermute(infoW))!=J(0,1,.)){
preDEN2=exp(quadcross(DUVEC,WVEC))
DEN2=DEN2+preDEN2
}


arssub1V ("alrank",DEN1,DEN2)
}


end






version 9
mata:
void arssub1V (transmorphic matrix input1a, real vector DEN1, real vector DEN2){ 
RVEC=st_data(.,input1a)



infoR=cvpermutesetup(RVEC) 
 
while ((RVEC=cvpermute(infoR))!=J(0,1,.)){
st_store(.,"__rorb",RVEC)
arssub2V ("__zorb",DEN1,DEN2)


}


}
end


version 9
mata:
void arssub2V (string scalar input2a,  real vector DEN1, real vector DEN2){ 


ZVEC=st_data(.,input2a)
RVEC=st_data(.,"__rorb")
//RVEC
//ZVEC

//RSUM=0

infoZ=cvpermutesetup(ZVEC) 
 
while ((ZVEC=cvpermute(infoZ))!=J(0,1,.)){
arssub3V(ZVEC,RVEC,DEN1,DEN2)
//ZVEC
//RVEC
}
}

end


version 9
mata:
void arssub3V (real vector zV, real vector rV, real vector DEN1, real vector DEN2){
RMUIK=st_data(1,"arunmuik")
//note where we are in loop; runmuik is the summed (ie, the actual) mu_ik
GUVEC=st_data(.,"_gunobs")
DUVEC=st_data(.,"_dunobs")



NUM1=exp(quadcross(GUVEC,zV))
NUM2=exp(quadcross(DUVEC,rV))

TESTST=quadcross(zV,rV)


MUIK=(TESTST)*((NUM1)/(DEN1))*((NUM2)/(DEN2))
//note that this is the summand of MU_ik; not MU_ik itself.



SIGMIK=((TESTST)-(RMUIK))^2*((NUM1)/(DEN1))*((NUM2)/(DEN2))
//note RMUIK is here the actual mu_ik, not the summand.

//(note case sensitivity...)
st_store(1,"muik",MUIK)
st_store(1,"sigmik",SIGMIK)



stata ("arssub4V muik sigmik")

}
end 


cap program drop arssub4V

program define arssub4V,rclass
version 9

syntax anything(name=numl)

gettoken one nlr: numl, parse(" ")
gettoken two nlr: nlr, parse(" ")




qui replace runmuik=runmuik + `one'
//note that this gives only one observation, which is what we want.


qui replace runsigmik=runsigmik + `two'


if (runmuik > usemui) | (usemui==. & runmuik !=0){
qui replace usemui=runmuik
qui replace usesigmi=runsigmik
}

if (runmui == usemui) & (runsigmik > usesigmi){
qui replace usemui=runmuik
qui replace usesigmi=runsigmik
}






end



//(dl:) cvpermute follows--omit, ultimately, and give code condt'l on Stata >= 10.
//! version 1.0.0  18aug2005
version 9.0

/*
	info = cvpermutesetup(V [, unique])

	real colvector cvpermute(info)


	Description of the contents of info:
		1.		1 -> unique
		2.		n = number of elements
		3.		m = current location
		4.		element 1
		5.		element 2
		...
		3+n.		element n

		3+n+1.		n-1
		3+n+2.		m = current location
		3+n+3.		element 2
		3+n+4.		element 3
		...
		etc.		until n-k = 2

	Let i0 be the base index - 1 for a group (i0=1 if the first 
	group, i0 = 3+n for the second, and so on), then 

		i0+1.		n = number of elements
		i0+2.		m = current location 
		i0+3.		element 1
		i0+4.		element 2
		...
		i0+3+n-1.	element n  and also next i0

	Example:

	Let V[] have N elements.

	N=2: (U,  2, #, x1, x2)
	N=3: (U,  3, #, x1, x2, x3,   2, #, y1, y2)
	N=4: (U,  4, #, x1, x2, x3, x4,   3, #, y1, y2, y3,   2, #, z1, z2)

	Ignoring U, 
		The 1st entry has N+2 elements
		The 2nd entry has N+2-1 elements
		The 3rd entry has N+2-2 element
		...
		The last entry has 4 elements

	total number of entries (ignoring U) is N+2 + N+2-1 + ... + 4.
*/


local UNIQUE 	1

local FIRST_n	2
local FIRST_m	3
local FIRST_el0	3

local FIRST_OFFSET	1

mata:

real rowvector cvpermutesetup(real colvector V, |real scalar unique)
{
	real scalar 	n, sum, i, isunique
	real rowvector	info
	real colvector	hold
	
	isunique = (args()==1 ? 0 : unique!=0)
	if ((n   = length(V))==0) return(J(1,3,0))
	if (n==1) {
		info = J(1, 4, 0)
		info[1] = isunique
		info[2] = info[3] = 1
		info[4] = V
		return(info)
	}

	/* 
		rows(info) = (4 + 5 + ... + n+2) + `FIST_OFFSET'

		The sum of (1,2,...k) is ((1+k)*k)/2
		Thus, the sum of (1, 2, ..., n+2) is ((3+n)*(n+2))/2
		and the from 4 is ((3+n)*(n+2))/2 - 6

		Alternative code:
		sum = 0 
		for (i=4; i<=n+2; i++) sum = sum + i 
	*/
	sum = round(((3+n)*(n+2))/2)-6
	info = J(1, sum+`FIRST_OFFSET', 0)

	info[`UNIQUE']  = isunique
	info[`FIRST_n'] = n
	info[`FIRST_m'] = 1

	if (isunique) info[|`FIRST_el0'+1\(`FIRST_el0'+n)|] = V'
	else {
		_sort(hold=V, 1)
		_transpose(hold)
		info[|`FIRST_el0'+1\(`FIRST_el0'+n)|] = hold
	}

	cvpermute_setupnext(info, `FIRST_OFFSET')
	return(info)
}


/* static */
void cvpermute_setupnext(
	real rowvector info,		// information vector
	real scalar i0)			// caller's base location in info[]
{
	real scalar n1, n2
	real scalar j0

	/*
	info :=               
                                               j0+1
                                                 |
                     i0+1   i0+3        i0+2+n1  |    j0+3
		       |      |            |     |      |
		(..., n1, m, x1, x2, ..., xn1,  n2, m, x1, x2, ..., xn2, ...)
                          |       |        |        |                |
                        i0+2    i0+4      j0      j0+2            j0+2+n2
	*/

	n1 = info[i0+1]
	if (n1==2) return

	j0 = i0+2+n1

	info[j0+1] = n2 = n1-1
	info[j0+2] = 1

	info[|j0+3\j0+2+n2|] = info[|i0+4\i0+2+n1|]

	cvpermute_setupnext(info, j0)
}

real colvector cvpermute(real rowvector info)
{
	real vector	res

	res = J(1, info[2], .)
	if (cvpermute_u(res, info, 1)) _transpose(res)
	else			     res = J(0, 1, .)
	return(res)
}


/* static */
real scalar cvpermute_u(real rowvector res, real rowvector info, real scalar i0)
{
	real scalar	i, n, m, mloc, lastel, firstel
	real scalar	nres
	real scalar	value

	nres = length(res)
	n = info[i0+1]
	m = info[mloc = i0+2]

	firstel = mloc + 1
	lastel  = mloc + n

	if (n==2) {
		if (m==1) {
			res[nres-1] = info[firstel]
			res[nres]   = info[lastel]
			info[mloc]  = 2
			return(1)
		}
		if (m==2) {
			if (!info[`UNIQUE']) {
				if (info[firstel]==info[lastel]) return(0)
			}
			res[nres-1] = info[lastel]
			res[nres]   = info[firstel]
			info[mloc]  = 3 
			return(1)
		}
		return(0)
	}

	if (n<=1) {		// n==1 arises only when there is 1 element
		if (n==0) return(0)
		if (m==1) {
			info[mloc] = 2 
			res[nres] = info[firstel]
			return(1)
		}
		return(0)
	}


	res[nres-n+1] = info[firstel]
	if (cvpermute_u(res, info, lastel)) return(1)

	value  = info[firstel]
	if (!info[`UNIQUE']) {
		for (i=mloc+m+1; i<=lastel; i++) {
			if (value != info[i]) {
				info[mloc] = i - mloc
				info[firstel] = info[i]
				info[i] = value
				res[nres-n+1] = info[firstel]
				cvpermute_setupnext(info, i0)
				return(cvpermute_u(res, info, lastel))
			}
		}
	}
	else {
		if ((i = mloc+m+1) <= lastel) {
			info[mloc] = i - mloc
			info[firstel] = info[i]
			info[i] = value
			res[nres-n+1] = info[firstel]
			cvpermute_setupnext(info, i0)
			return(cvpermute_u(res, info, lastel))
		}
	}
	return(0) 
}

end
